home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / specfn.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  14.2 KB  |  444 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module specfn)
  13.  
  14.    ;*********************************************************************
  15.    ;****************                                   ******************
  16.    ;**************** Macsyma Special Function Routines ******************
  17.    ;****************                                   ******************
  18.    ;*********************************************************************
  19.  
  20. (load-macsyma-macros rzmac)
  21. (load-macsyma-macros mhayat)
  22.  
  23. (defmacro mnumericalp (arg)
  24.       `(or (floatp ,arg) (and (or $numer $float) (integerp ,arg))))
  25.  
  26.         (comment Subtitle Polylogarithm Routines)
  27.  
  28. (declare-top(splitfile plylog))
  29.  
  30. (declare-top(special $zerobern ivars key-vars tlist)
  31.      (notype (li2numer flonum) (li3numer flonum)))
  32.  
  33. (defun lisimp (exp vestigial z)
  34.        vestigial ;; Ignored
  35.        (let ((s (simpcheck (car (subfunsubs exp)) z))
  36.          ($zerobern t)
  37.          (a))
  38.         (subargcheck exp 1 1 '$li)
  39.         (setq a (simpcheck (car (subfunargs exp)) z))
  40.         (or (cond ((zerop1 a) 0)
  41.               ((not (integerp s)) ())
  42.               ((= s 1) (m- `((%log) ,(m- 1 a))))
  43.               ((and (integerp a) (> s 1)
  44.                 (cond ((= a 1) ($zeta s))
  45.                   ((= a -1)
  46.                    (m*t (m1- `((rat) 1 ,(expt 2 (f- s 1))))
  47.                     ($zeta s))))))
  48.               ((= s 2) (li2simp a))
  49.               ((= s 3) (li3simp a)))
  50.         (eqtest (subfunmakes '$li (ncons s) (ncons a))
  51.             exp))))
  52.  
  53. (defun li2simp (arg)
  54.        (cond ((mnumericalp arg) (li2numer (float arg)))
  55.          ((alike1 arg '((rat) 1 2))
  56.           (m+t (m// ($zeta 2) 2)
  57.            (m*t '((rat simp) -1 2)
  58.             (m^ '((%log) 2) 2))))))
  59.  
  60. (defun li3simp (arg)
  61.        (cond ((mnumericalp arg) (li3numer (float arg)))
  62.          ((alike1 arg '((rat) 1 2))
  63.           (m+t (m* '((rat simp) 7 8) '(($zeta) 3))
  64.            (m*t (m// ($zeta 2) -2) (simplify '((%log) 2)))
  65.            (m*t '((rat simp) 1 6) (m^ '((%log) 2) 3))))))
  66.  
  67.  
  68. (declare-top (flonum x))
  69.  
  70. ;; Numerical evaluation for Chebyschev expansions of the first kind
  71.  
  72. (defun cheby (x chebarr)
  73.        (declare (flonum x))
  74.        (let ((Bn+2 0.0) (Bn+1 0.0))
  75.         (declare (flonum Bn+1 Bn+2))
  76.         (do ((i (fix (arraycall flonum chebarr 0)) (f1- i)))
  77.         ((< i 1) (-$ Bn+1 (*$ Bn+2 x)))
  78.         (setq Bn+2
  79.               (prog1 Bn+1 (setq Bn+1
  80.                     (+$ (arraycall flonum chebarr i)
  81.                         (-$ (*$ 2.0 x Bn+1)
  82.                         Bn+2))))))))
  83.  
  84. (defun cheby-prime (x chebarr)
  85.        (declare (flonum x))
  86.        (-$ (cheby x chebarr)
  87.        (*$ (arraycall flonum chebarr 1) .5)))
  88.  
  89. ;; These should really be calculated with minimax rational approximations.
  90. ;; Someone has done LI[2] already, and this should be updated; I haven't
  91. ;; seen any results for LI[3] yet.
  92.  
  93. (defun li2numer (x)
  94.        (cond ((= x 0.0) 0.0)
  95.          ((= x 1.0) 1.64493407)
  96.          ((= x -1.0) -.822467033)
  97.          ((< x -1.0) 
  98.           (-$ (+$ (chebyli2 (//$ x)) 1.64493407
  99.               (//$ (expt (log (-$ x)) 2) 2.0))))
  100.          ((not (> x .5)) (chebyli2 x))
  101.          ((< x 1.0)
  102.           (-$ 1.64493407 (*$ (log x) (log (-$ 1.0 x)))
  103.           (chebyli2 (-$ 1.0 x))))
  104.          (t (m+t (-$ 3.28986813 (//$ (expt (log x) 2) 2.0)
  105.              (li2numer (//$ x)))
  106.              (m*t (-$ (*$ 3.14159265 (log x))) '$%i)))))
  107.  
  108. (defun li3numer (x)
  109.        (cond ((= x 0.0) 0.0)
  110.          ((= x 1.0) 1.20205690)
  111.          ((< x -1.0)
  112.           (-$ (chebyli3 (//$ x)) (*$ 1.64493407 (log (-$ x)))
  113.           (//$ (expt (log (-$ x)) 3) 6.0)))
  114.          ((not (> x .5)) (chebyli3 x))
  115.          ((not (> x 2.0))
  116.           (let ((fac (*$ (expt (log x) 2) .5))) 
  117.            (m+t (+$ 1.20205690 
  118.                 (-$ (*$ (log x)
  119.                     (-$ 1.64493407 (chebyli2 (-$ 1.0 x))))
  120.                 (chebyS12 (-$ 1.0 x))
  121.                 (*$ fac
  122.                     (log (cond ((< x 1.0) (-$ 1.0 x))
  123.                            ((1-$ x)))))))
  124.             (cond ((< x 1.0) 0)
  125.                   ((m*t (*$ fac -3.14159265) '$%i))))))
  126.          (t (m+t (+$ (chebyli3 (//$ x)) (*$ 3.28986813 (log x))
  127.              (//$ (expt (log x) 3) -6.0))
  128.              (m*t (*$ -1.57079633 (expt (log x) 2)) '$%i)))))
  129.  
  130. ;(array *li2* flonum 15.)
  131. ;(array *li3* flonum 15.)
  132. ;(array *S12* flonum 18.)
  133.  
  134. (defvar *li2* (*array nil 'flonum 15.))
  135. (eval-when (load eval)
  136. (fillarray *li2* 
  137.        '(14.0       1.93506430 .166073033 2.48793229e-2 4.68636196e-3
  138.          1.0016275e-3 2.32002196e-4 5.68178227e-5 1.44963006e-5
  139.          3.81632946e-6 1.02990426e-6 2.83575385e-7 7.9387055e-8
  140.          2.2536705e-8 6.474338e-9))
  141. )
  142.  
  143. (defvar *li3* (*array nil 'flonum 15.))
  144. (eval-when (load eval)
  145. (fillarray *li3*
  146.        '(14.0       1.95841721 8.51881315e-2 8.55985222e-3 1.21177214e-3
  147.          2.07227685e-4 3.99695869e-5 8.38064066e-6 1.86848945e-6
  148.          4.36660867e-7 1.05917334e-7 2.6478920e-8 6.787e-9 
  149.          1.776536e-9 4.73417e-10))
  150. )
  151. (defvar *s12* (*array nil  'flonum 18.))
  152. (eval-when (load eval)
  153. (fillarray *S12*
  154.        '(17.0      1.90361778 .431311318 .100022507 2.44241560e-2
  155.          6.22512464e-3 1.64078831e-3 4.44079203e-4 1.22774942e-4
  156.          3.45398128e-5 9.85869565e-6 2.84856995e-6 8.31708473e-7
  157.          2.45039499e-7 7.2764962e-8 2.1758023e-8 6.546158e-9
  158.          1.980328e-9))
  159. )
  160. (defun chebyli2 (x)
  161.        (*$ x (cheby-prime (//$ (1+$ (*$ x 4.0)) 3.0) *li2*)))
  162.  
  163. ;              #+Maclisp '#,(get '*li2* 'array)
  164. ;              #+Franz (getd '*li2*)
  165. ;              #-(Or Maclisp Franz) (get-array-pointer '*li2*))))
  166.  
  167. (defun chebyli3 (x)
  168.        (*$ x (cheby-prime (//$ (1+$ (*$ 4.0 x)) 3.0) *li3*)))
  169. ;              #+Maclisp '#,(get '*li3* 'array)
  170. ;              #+Franz (getd '*li3*)
  171. ;              #-(or Maclisp Franz) (get-array-pointer '*li3*))))
  172.  
  173. (defun chebyS12 (x)
  174.   (*$ (//$ (expt x 2) 4.0)
  175.       (cheby-prime (//$ (1+$ (*$ 4.0 x)) 3.0) *s12*)))
  176. ;           #+Maclisp '#,(get '*S12* 'array)
  177. ;           #+Franz (getd '*S12*)
  178. ;           #-(or Maclisp Franz) (get-array-pointer '*S12*))))
  179.  
  180. (declare-top(notype x))
  181.  
  182.         (comment Subtitle Polygamma Routines)
  183.  
  184. (declare-top(splitfile plygam))
  185.  
  186. ;; gross efficiency hack, exp is a function of *k*, *k* should be mbind'ed
  187.  
  188. (defun msum (exp lo hi)
  189.   (if (< hi lo)
  190.       0
  191.       (let ((sum 0))
  192.     (do ((*k* lo (f1+ *k*)))
  193.         ((> *k* hi) sum)
  194.        (declare (special *k*))
  195.       (setq sum (add2 sum (meval exp)))))))
  196.  
  197.  
  198. (defun pole-err (exp)
  199.        (declare (special errorsw))
  200.        (cond (errorsw (throw 'errorsw t))
  201.          (t (merror "Pole encountered in: ~M" exp)
  202.         )))
  203.  
  204.  
  205. (declare-top (special
  206.       $maxpsiposint $maxpsinegint $maxpsifracnum $maxpsifracdenom)
  207.      (fixnum 
  208.       $maxpsiposint $maxpsinegint $maxpsifracnum $maxpsifracdenom)
  209.      (*lexpr $diff))
  210.  
  211. (defprop $psi psisimp specsimp)
  212.  
  213. (mapcar (function (lambda (var val)
  214.               (and (not (boundp var)) (set var val))))
  215.     '($maxpsiposint $maxpsinegint $maxpsifracnum $maxpsifracdenom)
  216.     '(20. -10. 4 4))
  217.  
  218. (defun psisimp (exp a z)
  219.        (let ((s (simpcheck (car (subfunsubs exp)) z)))
  220.         (subargcheck exp 1 1 '$psi)
  221.         (setq a (simpcheck (car (subfunargs exp)) z))
  222.         (and (integerp a) (< a 1) (pole-err exp))
  223.         (eqtest (psisimp1 s a) exp)))
  224.  
  225. ;; This gets pretty hairy now.
  226.  
  227. (defun psisimp1 (s a)
  228.  (let ((*k*))
  229.       (declare (special *k*))
  230.       (or
  231.        (and (not $numer) (not $float) (integerp s) (> s -1)
  232.     (cond
  233.      ((integerp a)
  234.       (and (not (> a $maxpsiposint))    ; integer values
  235.            (m*t (expt -1 s) (factorial s)
  236.             (m- (msum (inv (m^t '*k* (f1+ s))) 1 (f1- a))
  237.             (cond ((zerop s) '$%gamma)
  238.                   (($zeta (f1+ s))))))))
  239.      ((or (not (ratnump a)) (ratgreaterp a $maxpsiposint)) ())
  240.      ((ratgreaterp a 0)
  241.       (cond
  242.        ((ratgreaterp a 1)
  243.         (let* ((int ($entier a))    ; reduction to fractional values
  244.            (frac (m-t a int)))
  245.           (m+t
  246.            (psisimp1 s frac)
  247.            (if (> int $maxpsiposint)
  248.                (subfunmakes '$psi (ncons s) (ncons int))
  249.                (m*t (expt -1 s) (factorial s)
  250.                 (msum (m^t (m+t (m-t a int) '*k*)
  251.                        (f1- (f- s)))
  252.                   0 (f1- int)))))))
  253.        ((= s 0)
  254.         (let ((p (cadr a)) (q (caddr a)))
  255.          (cond
  256.           ((or (greaterp p $maxpsifracnum)
  257.            (greaterp q $maxpsifracdenom) (bigp p) (bigp q)) ())
  258.           ((and (= p 1)
  259.         (cond ((= q 2)
  260.                (m+ (m* -2 '((%log) 2)) (m- '$%gamma)))
  261.               ((= q 3)                            
  262.                (m+ (m* '((rat simp) -1 2)
  263.                    (m^t 3 '((rat simp) -1 2)) '$%pi)
  264.                (m* '((rat simp) -3 2) '((%log) 3))
  265.                (m- '$%gamma)))
  266.               ((= q 4)
  267.                (m+ (m* '((rat simp) -1 2) '$%pi)
  268.                (m* -3 '((%log) 2)) (m- '$%gamma))))))
  269.           ((and (= p 2) (= q 3))
  270.            (m+ (m* '((rat simp) 1 2)
  271.                (m^t 3 '((rat simp) -1 2)) '$%pi)
  272.            (m* '((rat simp) -3 2) '((%log) 3))
  273.            (m- '$%gamma)))
  274.           ((and (= p 3) (= q 4))
  275.            (m+ (m* '((rat simp) 1 2) '$%pi)
  276.            (m* -3 '((%log) 2)) (m- '$%gamma)))
  277.           ;; Gauss's Formula
  278.           ((let ((f (m* `((%cos) ,(m* 2 a '$%pi '*k*))
  279.                 `((%log) ,(m-t 2 (m* 2 `((%cos) 
  280.                              ,(m//t (m* 2 '$%pi '*k*)
  281.                                 q))))))))
  282.             (m+t (msum f 1 (f1- (// q 2)))
  283.              (let ((*k* (// q 2)))
  284.                (declare (special *k*))
  285.                   (m*t (meval f)
  286.                    (cond ((oddp q) 1)
  287.                      ('((rat simp) 1 2)))))
  288.              (m-t (m+ (m* '$%pi '((rat simp) 1 2)
  289.                       `((%cot) ((mtimes simp) ,a $%pi)))
  290.                   `((%log) ,q)
  291.                   '$%gamma))))))))
  292.        ((alike1 a '((rat) 1 2))
  293.         (m*t (expt -1 (f1+ s)) (factorial s)
  294.          (f1- (expt 2 (f1+ s))) (simplify ($zeta (f1+ s)))))))
  295.      ((ratgreaterp a $maxpsinegint)  ;;; Reflection Formula
  296.       (m*t
  297.        (^ -1 s)
  298.        (m+t (m+t (psisimp1 s (m- a))
  299.              (let ((dif ($diff `((%cot) ,(m* '$%pi '$z)) '$z s))
  300.                ($z (m-t a)))
  301.              (declare (special $z))
  302.              (meval dif)))
  303.         (m*t (factorial s) (m^t (m-t a) (f1- (f- s)))))))))
  304.        (subfunmakes '$psi (ncons s) (ncons a)))))
  305.  
  306.  
  307.         (comment Subtitle Polygamma Tayloring Routines)
  308.  
  309. ;; These routines are specially coded to be as fast as possible given the
  310. ;; current $TAYLOR; too bad they have to be so ugly.
  311.  
  312. (declare-top(special var subl last sign last-exp))
  313.  
  314. (defun expgam-fun (pw temp)
  315.        (setq temp (get-datum (get-key-var (car var))))
  316.        (let-pw temp pw
  317.            (pstimes
  318.         (let-pw temp (e1+ pw)
  319.             (psexpt-fn (getexp-fun '(($PSI) -1) var (e1+ pw))))
  320.         (make-ps var (ncons pw) '(((-1 . 1) 1 . 1))))))
  321.  
  322. (defun expplygam-funs (pw subl l) ; l is a irrelevant here
  323.        (setq subl (car subl))
  324.        (if (or (not (integerp subl)) (lessp subl -1))
  325.        (tay-err "Unable to expand at a subscript in")
  326.        (prog ((e 0) (sign 0) npw)
  327.         (declare (flonum npw) (fixnum e) #-Multics (fixnum sign))
  328.         (setq npw (//$ (float (car pw)) (float (cdr pw))))
  329.         (setq
  330.          l (cond ((= subl -1)
  331.               `(((1 . 1) . ,(prep1 '((mtimes) -1 $%gamma)))))
  332.              ((= subl 0)
  333.               (cons '((-1 . 1) -1 . 1)
  334.                 (if (> 0.0 npw) ()
  335.                 `(((0 . 1)
  336.                    . ,(prep1 '((mtimes) -1 $%gamma)))))))
  337.              (T (setq last (factorial subl))
  338.             `(((,(f- (f1+ subl)) . 1)
  339.                ,(times (^ -1 (f1+ subl))
  340.                    (factorial subl)) . 1))))
  341.          e (if (< subl 1) (f- subl) -1)
  342.          sign (if (< subl 1) -1 (^ -1 subl)))
  343.       a (setq e (f1+ e) sign (f- sign))
  344.         (if (greaterp e npw) (return l)
  345.         (rplacd (last l)
  346.             `(((,e . 1)
  347.                . ,(rctimes (rcplygam e)
  348.                        (prep1 ($zeta (f+ (f1+ subl) e))))))))
  349.         (go a))))
  350.  
  351. (defun rcplygam (k)
  352.        (declare (fixnum k) )
  353.        (cond ((= subl -1) (cons sign k))
  354.          ((= subl 0) (cons sign 1))
  355.          (t (prog1 (cons (times sign last) 1)
  356.           
  357.                (setq last
  358.                  (*quo (times last (plus subl (add1 k)))
  359.                    (add1 k)))))))
  360.  
  361. (defun plygam-ord (subl)
  362.        (if (equal (car subl) -1) (ncons (rcone))
  363.        `((,(f- (f1+ (car subl))) . 1))))
  364.  
  365. (defun plygam-pole (a c func)
  366.        (if (rcmintegerp c)
  367.        (let ((ps (get-lexp (m- a (rcdisrep c)) () t)))
  368.         (rplacd (cddr ps) (cons `((0 . 1) . ,c) (cdddr ps)))
  369.         (if (atom func) (gam-const a ps func)
  370.             (plygam-const a ps func)))
  371.        (prep1 (simplifya
  372.            (if (atom func) `((%GAMMA) ,(rcdisrep c))
  373.                `((MQAPPLY) ,func ,(rcdisrep c)))
  374.            () ))))
  375.  
  376. (defun gam-const (a arg func)
  377.    (let ((const (ps-lc* arg)) (arg-c))
  378.     (ifn (rcintegerp const)
  379.          (taylor2 (diff-expand `((%GAMMA) ,a) tlist))
  380.          (setq const (car const))
  381.          ;; Try to get the datum
  382.          (if (pscoefp arg)
  383.          (setq arg-c (get-lexp (m+t a (minus const)) (rcone)
  384.                        (signp le const))))
  385.          (if (and arg-c (not (psp arg-c)))    ; must be zero
  386.          (taylor2 (simplify `((%GAMMA) ,const)))
  387.          (let ((datum (get-datum (get-key-var
  388.                       (gvar (or arg-c arg)))))
  389.                (ord (if arg-c (le (terms arg-c))
  390.                 (le (n-term (terms arg))))))
  391.             (setq func (current-trunc datum))
  392.             (if (greaterp const 0)
  393.             (pstimes 
  394.              (let-pw datum (e- func ord)
  395.                  (expand (m+t a (minus const)) '%GAMMA))
  396.              (let-pw datum (e+ func ord)
  397.                  (tsprsum (m+t a (m-t '%%taylor-index%%))
  398.                       `(%%taylor-index%% 1 ,const)
  399.                       '%PRODUCT)))
  400.             (pstimes 
  401.              (expand (m+t a (minus const)) '%GAMMA)
  402.              (let-pw datum (e+ func ord)
  403.                  (psexpt 
  404.                   (tsprsum (m+t a '%%taylor-index%%)
  405.                        `(%%taylor-index%% 0
  406.                            ,(minus (add1 const))) '%PRODUCT)
  407.                   (rcmone))))))))))
  408.  
  409. (defun plygam-const (a arg func)
  410.    (let ((const (ps-lc* arg)) (sub (cadr func)))
  411.     (cond 
  412.      ((or (not (integerp sub)) (< sub -1))
  413.       (tay-err "Unable to expand at a subscript in"))
  414.      ((not (rcintegerp const))
  415.       (taylor2 (diff-expand `((MQAPPLY) ,func ,a) tlist)))
  416.      (T (setq const (car const))
  417.         (psplus
  418.          (expand (m+t a (f- const)) func)
  419.          (if (> const 0)
  420.          (pstimes
  421.           (cons (times (^ -1 sub) (factorial sub)) 1)
  422.           (tsprsum `((MEXPT) ,(m+t a (m-t '%%taylor-index%%)) ,(f- (f1+ sub)))
  423.                `(%%taylor-index%% 1 ,const) '%SUM))
  424.          (pstimes
  425.           (cons (times (^ -1 (f1+ sub)) (factorial sub)) 1)
  426.           (tsprsum `((MEXPT) ,(m+t a '%%taylor-index%%) ,(f- (f1+ sub)))
  427.                `(%%taylor-index%% 0 ,(f- (f1+ const))) '%SUM))))))))
  428.  
  429. (declare-top(unspecial var subl last sign last-exp))
  430.  
  431. ;; Not done correctly
  432. ;;
  433. ;; (defun beta-trans (argl funname)
  434. ;;   funname ;ignored
  435. ;;   (let ((sum (m+t (car argl) (cadr argl))) (PSI[-1] '(($PSI ARRAY) -1)))
  436. ;;     (if (zerop sum) (unfam-sing-err)
  437. ;;     (taylor2 `((MTIMES)
  438. ;;            ((MEXPT) $%E ((MPLUS) ((MQAPPLY) ,PSI[-1] ,(car argl))
  439. ;;                      ((MQAPPLY) ,PSI[-1] ,(cadr argl))
  440. ;;                      ((MTIMES) -1
  441. ;;                            ((MQAPPLY) ,PSI[-1] ,sum))))
  442. ;;            ((MPLUS) ((MEXPT) ,(car argl) -1)
  443. ;;                 ((MEXPT) ,(cadr argl) -1)))))))
  444.